Impact of water deficit on growth, productivity, and water use efficiency in potato genotypes (Solanum tuberosum L.)

1 Setup

setup <- tempfile(fileext = ".r")
readLines(con = 'https://inkaverse.com/setup.r') |> 
  gsub("echo = FALSE", "echo = TRUE", x = _) |>
  writeLines(setup)
source(setup)

library(gtsummary)
library(factoextra)
library(corrplot)

2 Googlesheets connect

url <- "https://docs.google.com/spreadsheets/d/1dfgpmCKdPmxRHozrZp0iE_xMGsvKTIcztDpMWYSEGaY"
gs <- as_sheets_id(url)
# browseURL(url)

3 Trait rename

geno <- c("CIP720088" = "G01"
          , "CIP392797.22" = "G02"
          , "CIP397077.16" = "G03"
          , "CIP398192.213" = "G04"
          , "CIP398180.612" = "G05"
          , "CIP398208.704" = "G06"
          , "CIP398098.119" = "G07"
          , "CIP398190.89" = "G08"
          , "CIP398192.592" = "G09"
          , "CIP398201.510" = "G10"
          , "CIP398203.244" = "G11"
          , "CIP398203.5" = "G12"
          , "CIP398208.219" = "G13"
          , "CIP398208.33" = "G14"
          , "CIP398208.620" = "G15")

traits <- c("SPAD_29" = "Chlorophyll concentration (SPAD) at 29 dap"
            , "SPAD_59" = "Chlorophyll concentration (SPAD) at 59 dap"
            , "SPAD_76" = "Chlorophyll concentration (SPAD) at 76 dap"
            , "SPAD_83" = "Chlorophyll concentration (SPAD) at 83 dap"
            , "HGT" = "Plant height (cm)"
            , "RWC" = "Relative water content (%)"
            , "LOP" = "Leaf osmotic potential (MPa)"
            , "LDW" = "Leaf dry weight (g)"
            , "SDW" = "Stem dry weight (g)"
            , "RDW" = "Root dry weight (g)"
            , "TDW" = "Tuber dry weight (g)"
            , "NTUB" = "Tuber number (N°)" 
            , "TRS" = "Total transpiration (mL)"
            , "LFA" = "Leaf area (cm^2^)"
            , "RTL" = "Root length (cm)"
            , "TDB" = "Total dry biomass (g)"
            , "HI"  = "Harvest index (HI)" 
            , "SLA" = "Specific leaf area (cm^2^/g)"
            , "RCC" = "Relative chlorophyll content (RCC)"
            , "WUE_B" = "Biomass water use efficiency (g/L)"
            , "WUE_T" = "Tuber water use efficiency (g/L)"
            )

4 Import data

fb <- gs %>% 
  range_read("fb") %>% 
  rename_with("tolower") %>%
  rename_with(~gsub("\\s+|\\.", "_", .)) %>% 
  mutate(treat = case_when(treat %in% "wellwater" ~ "WW"
                           , treat %in% "drydown" ~ "WD")) %>% 
  select(block, treat, genotype,
         spad_29 = spad_29dap, 
         spad_59 = spad_59dap, 
         spad_76 = spad_76dap, 
         spad_83 = spad_83dap,
         hgt = hgt_86dap,
         rwc = rwc_84dap,
         lop = op_84dap,
         ldw = leafdw,
         sdw = stemdw,
         rdw = rootdw,
         tdw = tubdw,
         ntub,
         trs = ttrns,
         lfa = la,
         rtl = rlg
         ) %>% 
  mutate(tdb = (ldw+sdw+rdw+tdw),
         hi = tdw/(ldw+sdw+rdw+tdw),
         sla = lfa/ldw,
         rcc = (spad_83/lfa),
         wue_b = (ldw+sdw+rdw+tdw)/trs,
         wue_t = tdw/trs
         ) %>% 
  mutate(gnt = recode(genotype, !!!geno), .after = genotype) %>% 
  select(genotype, gnt, treat, block, everything()) %>% 
  mutate(across(where(is.character), as.factor)) 

# str(fb)

fb %>% web_table()

5 Abbreviations

gs %>% 
  range_read("var") %>% 
  filter(include == "x") %>%
  select(Variable, Abbreviation) %>%
  kable()
Abbreviations
Variable Abbreviation
Soil Plant Analysis Development SPAD
Plant height HGT
Relative water content RWC
Leaf osmotic potential LOP
Leaf dry weight LDW
Stem dry weight SDW
Root dry weight RDW
Tuber dry weight TDW
Tuber number NTUB
Root length RTL
Total transpiration TRS
Leaf area LFA
Total dry biomass TDB
Harvest index HI
Specific leaf area SLA
Water use efficiency WUE
Tuber water use efficiency TWUE

6 Table 1

tab <- gs %>% 
  range_read("gnt") %>% 
  dplyr::select(Number, Genotypes, Adaptability, "Growing period" = GPL, "Heat tolerance" = Heat, "Dry matter (%)")
  
# tab %>% sheet_write(gs, data = .,  sheet = "tab1")

tab %>% web_table()

7 Table 2

tab <- fb %>%
  select(!c(block, gnt, genotype)) %>% 
  tbl_summary(
    by = treat, 
    statistic = list(all_continuous() ~ "{mean} ± {sd}"),
    missing = "no") %>% 
  add_p() %>% 
  bold_labels() %>% 
  as_tibble() %>% 
  mutate(`**Characteristic**` = recode(
    `**Characteristic**`
    , "__spad_29__" = "Chlorophyll concentration (SPAD) at 29 dap"
    , "__spad_59__" = "Chlorophyll concentration (SPAD) at 59 dap"
    , "__spad_76__" = "Chlorophyll concentration (SPAD) at 76 dap"
    , "__spad_83__" = "Chlorophyll concentration (SPAD) at 83 dap"
    , "__hgt__" = "Plant height (cm)"
    , "__rwc__" = "Relative water content (%)"
    , "__lop__" = "Leaf osmotic potential (MPa)"
    , "__ldw__" = "Leaf dry weight (g)"
    , "__sdw__" = "Stem dry weight (g)"
    , "__rdw__" = "Root dry weight (g)"
    , "__tdw__" = "Tuber dry weight (g)"
    , "__ntub__" = "Tuber number (N°)" 
    , "__trs__" = "Total transpiration (mL)"
    , "__lfa__" = "Leaf area (cm^2^)"
    , "__rtl__" = "Root length (cm)"
    , "__tdb__" = "Total dry biomass (g)"
    , "__hi__"  = "Harvest index (HI)" 
    , "__sla__" = "Specific leaf area (cm^2^/g)"
    , "__rcc__" = "Relative chlorophyll content (RCC)"
    , "__wue_b__" = "Biomass water use efficiency (g/L)"
    , "__wue_t__" = "Tuber water use efficiency (g/L)"
    )) %>% 
  rename(Variable = `**Characteristic**`
         , `Water deficit` = `**WD**, N = 75`
         , `Well-Watered`= `**WW**, N = 75`
         , `p-value` = `**p-value**`) 

# tab %>% sheet_write(gs, data = .,  sheet = "tab2")

tab %>% web_table()

8 Figure 1

fts <- gs %>% 
  range_read("ftsw") %>%
  filter(Treatment != "preharvest") %>%
  tidyr::gather(key = day, value = fts, -ID, -Genotype, -Treatment)

model <- fts %>% 
  inti::mean_comparison(response = "fts"
                        , model_factors = "Treatment*day"
                        , c("Treatment", "day"))

plt1 <- model$comparison %>% 
  mutate(Treatment = case_when(
    Treatment == "drydown" ~ "Water Deficit (WD)"
    , Treatment == "wellwater" ~ "Well Watered (WW)"
  )) %>% 
  mutate(across(day, as.numeric)) %>% 
  plot_smr(
  type = "line", color = F
  , x = "day"
  , y = "fts"
  , group = "Treatment"
  , ylab = "Fraction of transpirable soil water"
  , xlab =  "Days"
  , glab = "Treatment"
  , legend = "none"
  , ylimits = c(0, 1.08, 0.1)
  , error =  "ste"
  ) +
  scale_x_continuous(breaks = c(0:100)*2
                     , limits = c(34, 88))

trns <- gs %>% 
  range_read("trans") %>%
  filter(Treatment != "preharvest") %>%
  tidyr::gather(key = day, value = trans, -ID, -Genotype, -Treatment) %>%
  filter(day != "TOTAL") %>%
  drop_na()


model <- trns %>% 
  inti::mean_comparison(response = "trans"
                        , model_factors = "Treatment*day"
                        , c("Treatment", "day"))

plt2 <- model$comparison %>% 
    mutate(Treatment = case_when(
    Treatment == "drydown" ~ "Water Deficit (WD)"
    , Treatment == "wellwater" ~ "Well Watered (WW)"
  )) %>% 
  mutate(across(day, as.numeric)) %>% 
  plot_smr(type = "line", color = F
            , x = "day"
            , y = "trans"
            , group = "Treatment"
            , ylab = "Transpiration (mL)"
            , xlab =  "Days"
            , glab = "Irrigation"
            , ylimits = c(0, 700, 100)
            , error = "ste"
            ) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

plot <- ggdraw(xlim = c(0, 0.5), ylim = c(0,0.9)) +
  draw_plot(plt2, width = 0.5, height = 0.43 , x = 0.0, y = 0.47 ) +
  draw_plot(plt1, width = 0.495, height = 0.5, x = 0.005, y = 0.0 ) +
          draw_plot_label(
            label = c("a", "b"),
            x = c(0.465, 0.465),
            y = c(0.79, 0.49))

ggsave2(plot = plot, "files/Fig1.tiff", width = 18, height = 15, units = "cm")

"files/Fig1.jpg" %>% include_graphics()

9 Figure 2

# tdw

model <- fb %>% 
  inti::mean_comparison(response = "tdw"
                        , model_factors = "block + gnt*treat"
                        , c("gnt","treat"))

tdw <- model$comparison %>% 
  plot_smr(type = "bar"
          , color = F
          , x = "gnt"
          , xlab = ""
          , y = "tdw"
          , ylab = "Tuber dry weight (g)"
          , group = "treat"
          , glab = "Treatment"
          , legend = "none"
          , error = "ste"
          , ylimits = c(0, 100, 20)
          )

# rcc

model <- fb %>% 
  inti::mean_comparison(response = "rcc"
                        , model_factors = "block + gnt*treat"
                        , c("gnt","treat"))

rcc <- model$comparison %>% 
  plot_smr(type = "bar"
        , color = F
        , x = "gnt"
        , xlab =  ""
        , y = "rcc"
        , ylab = "Relative chlorophyll content"
        , ylimits = c(0, 0.1, 0.02)
        , group = "treat"
        , glab = "Treatment"
        , legend = "none"
        , error = "ste"
        ) 

# hi

model <- fb %>% 
  inti::mean_comparison(response = "hi"
                        , model_factors = "block + gnt*treat"
                        , c("gnt","treat"))

hi <- model$comparison %>% 
  plot_smr(type = "bar"
        , color = F
        , x = "gnt"
        , xlab =  ""
        , y = "hi"
        , ylab = "Harvest Index"
        , group = "treat"
        , glab = "Treatment"
        , legend = "none"
        , ylimits = c(0, 1, 0.2)
        , error = "ste"
        ) 

# wue_t

model <- fb %>% 
  inti::mean_comparison(response = "wue_t"
                        , model_factors = "block + gnt*treat"
                        , c("gnt","treat"))

wue_t <- model$comparison %>% 
  plot_smr(type = "bar", color = F
        , x = "gnt"
        , xlab =  ""
        , y = "wue_t"
        , ylab = "Tuber water use efficiency (g/L)"
        , group = "treat"
        , glab = "Treatment"
        , legend = "none"
        , ylimits = c(0, 10, 2)
        , error = "ste"
        ) 

# spad_29

model <- fb %>% 
  inti::mean_comparison(response = "spad_29"
                        , model_factors = "block + gnt*treat"
                        , c("gnt","treat"))

spad_29 <- model$comparison %>% 
    mutate(treat  = case_when(
    treat  == "WD" ~ "Water Deficit (WD)"
    , treat  == "WW" ~ "Well Watered (WW)"
    )) %>%
  plot_smr(type = "line"
        , x = "gnt"
        , xlab =  ""
        , y = "spad_29"
        , ylab = "SPAD at 29 dap"
        , group = "treat"
        , glab = "Treatment"
        , legend = "top"
        , ylimits = c(30, 70, 5)
        , error = "ste"
        , color = F
        )

# spad_83

model <- fb %>% 
  inti::mean_comparison(response = "spad_83"
                        , model_factors = "block + gnt*treat"
                        , c("gnt","treat"))

spad_83 <- model$comparison %>% 
  plot_smr(type = "line", color = F
        , x = "gnt"
        , xlab =  ""
        , y = "spad_83"
        , ylab = "SPAD at 83 dap"
        , group = "treat"
        , glab = "Treatment"
        , legend= "none"
        , ylimits = c(30, 70, 5)
        , error = "ste"
        ) 

plot <- ggdraw(xlim = c(0.0, 1), ylim = c(0.0, 1)) +
  draw_plot(spad_29 , width = 0.49, height = 0.35, x = 0.0, y = 0.6) +
  draw_plot(spad_83, width = 0.49, height = 0.3, x = 0.5, y = 0.6) +
  draw_plot(rcc,  width = 0.49, height = 0.3, x = 0.0, y = 0.3) +
  draw_plot(tdw, width = 0.49, height = 0.3, x = 0.5, y = 0.3) +
  draw_plot(hi , width = 0.49, height = 0.3, x = 0.0, y = 0.0) +
  draw_plot(wue_t, width = 0.49, height = 0.3, x = 0.5, y = 0.0) +
          draw_plot_label(
            label = c("a", "b", "c", "d", "e", "f"),
            x = c(0.06, 0.56, 0.06, 0.56, 0.06, 0.56),
            y = c(0.89, 0.89, 0.59, 0.59, 0.29, 0.29))

ggsave(plot = plot, "files/Fig2.tiff", units = "cm", width = 28, height = 25)

"files/Fig2.jpg" %>% include_graphics()

10 Figure 3

cr <- fb %>% 
  group_by(gnt, treat) %>% 
  summarise_all(list(mean), na.rm = TRUE) %>% 
  as.data.frame() %>% 
  select(!c( genotype, block)) %>% 
  unite(gnt, treat, col = "rn") %>% 
  column_to_rownames("rn") %>% 
  rename_with(toupper, where(is.numeric)) 

cp <- heatmaply::heatmaply_cor(
  cor(cr),
  k_col = 5, 
  dendrogram = c("column"),  
  grid_gap = 1,
  cellnote = round(cor(cr),2),
  cellnote_textposition = "middle center",
  cellnote_size = 10,
  file = "files/correlation.pdf"
  )

fig <- magick::image_read_pdf("files/correlation.pdf") %>% 
  grid::rasterGrob() %>%
  ggsave("files/Fig3.tiff", plot =  ., width = 8.39, height = 5.26)
  
unlink("files/correlation.pdf")

"files/Fig3.jpg" %>% include_graphics()

Where: Chlorophyll Concentration (SPAD), Plant height (HGT; cm), Relative water content (RWC; %), Leaf osmotic potential (LOP; MPa), Leaf dry weight (LDW; g), Stem dry weight (SDW; g), Root dry weight (RDW; g), Tuber dry weight (TDW; g), Tuber number (NTUB; N°), Total transpiration (TRS; mL), Leaf area (LFA; cm2), Root length (RTL; cm), Total dry biomass (TDB; g), Harvest index (HI), Specific leaf area (SLA; cm2g-1), Relative chlorophyll content (RCC), Biomass water use efficiency (WUEB; gL-1), Tuber water use efficiency (WUET; gL-1).

11 Figure 4

mvd <- fb %>%
  select(!c(block, gnt)) %>%
  group_by(treat, genotype) %>%
  summarise(across(everything(), ~mean(.x, na.rm = T))) %>% 
  mutate(coln = paste(genotype, treat,  sep = "_")) %>%
  column_to_rownames("coln") %>%
  select(!genotype) %>% 
  rename_with(toupper, where(is.numeric))

pca <- PCA(mvd, graph = F, scale.unit = TRUE, quali.sup = 1)

# PCA

tmp <- tempfile(fileext = ".png")
ppi <- 300
png(tmp, width=8*ppi, height=8*ppi, res=ppi)
plot.PCA(pca,
         choix="var",
         title="",
         autoLab = "auto", 
         shadowtext = T)
graphics.off()

pt1 <- png::readPNG(tmp) %>%
  grid::rasterGrob()

# sink("files/pca.txt")
# # Results
# summary(pca, nbelements = Inf, nb.dec = 2)
# # Correlation de dimensions
# dimdesc(pca)
# sink()

# Analysis de Hierarchical Clustering

clus <- HCPC(pca, nb.clust=-1, graph = F)

tmp <- tempfile(fileext = ".png")
ppi <- 300
png(tmp, width=8*ppi, height=8*ppi, res=ppi)
plot.HCPC(clus,
          choice = "map",
          legend = list(bty = "y", x = "topright"), 
          draw.tree = F)
graphics.off()

pt2 <- png::readPNG(tmp) %>%
  grid::rasterGrob()

# sink("files/clu.txt")
# clus$call$t$tree
# clus$desc.ind
# clus$desc.var
# sink()

plot <- plot_grid(pt1, pt2, nrow = 1, labels = "auto")

ggsave2(plot = plot, "files/Fig4.tiff", units = "cm", width = 20, height = 11)

"files/Fig4.jpg" %>% include_graphics()

Where: Chlorophyll Concentration (SPAD), Plant height (HGT; cm), Relative water content (RWC; %), Leaf osmotic potential (LOP; MPa), Leaf dry weight (LDW; g), Stem dry weight (SDW; g), Root dry weight (RDW; g), Tuber dry weight (TDW; g), Tuber number (NTUB; N°), Total transpiration (TRS; mL), Leaf area (LFA; cm2), Root length (RTL; cm), Total dry biomass (TDB; g), Harvest index (HI), Specific leaf area (SLA; cm2g-1), Relative chlorophyll content (RCC), Biomass water use efficiency (WUEB; gL-1), Tuber water use efficiency (WUET; gL-1).

12 Supplementary Table 1

vars <- names(fb)

model <- 5:length(vars) %>% map(function(x) {
  
    fb %>% 
      H2cal(trait = vars[x]
            , gen.name = "genotype"
            , rep.n = 4
            , fixed.model = "0 + (1|block) + genotype"
            , random.model = "1 + (1|block) + (1|genotype)"
            , emmeans = T
            , plot_diag = T
            , outliers.rm = T
            )
})


tabsmr <- 1:length(model) %>% map(function(x) {
      model[[x]][["tabsmr"]] 
      }) %>% 
  bind_rows() %>% 
  mutate(trait = recode(trait, !!!traits)) %>% 
  select(trait, mean, std, min, max, V.g, V.e, repeatability) %>% 
  mutate(across(where(is.numeric), ~ round(.x, digits = 2))) %>% 
  mutate(trait = toupper(trait)) %>% 
  mutate(trait = recode(trait, !!!traits)) %>% 
  rename(Trait = trait)
  
tabsmr  %>% web_table()

# tabsmr %>% sheet_write(gs, data = .,  sheet = "tabS1")

13 Supplementary Figure 1

var <- get_pca_var(pca)

tmp <- tempfile(fileext = ".png")
ppi <- 300
png(tmp, width=3.8*ppi, height=10*ppi, res=ppi)
corrplot(var$cor, 
         method="number",
         tl.col="black", 
         tl.srt=45,)
graphics.off()

pt1 <- png::readPNG(tmp) %>%
  grid::rasterGrob(interpolate = TRUE)

pt2 <- fviz_eig(pca, 
                addlabels=TRUE,
                hjust = 0.05,
                barfill="white",
                barcolor ="darkblue",
                linecolor ="red") + 
  ylim(0, 50) + 
  labs(
    title = "PCA - percentage of explained variances",
    y = "Variance (%)") +
  theme_minimal()

pt3 <- fviz_contrib(pca,
                     choice = "var", 
                     axes = 1, 
                     top = 10,
                     fill="white",
                     color ="darkblue",
                     sort.val = "desc") +
  ylim(0, 18) + 
  labs(title = "Dim 1 - variables contribution") 

pt4 <- fviz_contrib(pca,
                     choice = "var", 
                     axes = 2, 
                     top = 10,
                     fill="white",
                     color ="darkblue",
                     sort.val = "desc") +
  ylim(0, 18) + 
  labs(title = "Dim 2 - variables contribution") 

plot <- ggdraw(xlim = c(0.0, 1.0), ylim = c(0, 1.0))+
  draw_plot(pt1,  width = 0.4, height = 0.99, x = 0.6, y = 0.0) +  
  draw_plot(pt2,  width = 0.6, height = 0.34, x = 0.03, y = 0.66) +
  draw_plot(pt3, width = 0.6, height = 0.34, x = 0.03, y = 0.33) + 
  draw_plot(pt4, width = 0.6, height = 0.34, x = 0.03, y = 0.0) +
        draw_plot_label(
    label = c("a", "b", "c", "d"),
    x = c(0.005, 0.005, 0.005, 0.65),
    y = c(0.999, 0.67, 0.34, 0.999))

ggsave2(plot = plot, "files/SFig1.tiff", height = 25, width = 30, units = "cm")

"files/SFig1.jpg" %>% include_graphics()

14 Supplementary Figure 2

fig <- magick::image_read("files/potatos.png") %>% 
  grid::rasterGrob() %>%
  ggsave("files/SFig2.tiff", plot =  ., width = 15, height = 15)

"files/SFig2.jpg" %>% include_graphics()

15 Convert figures to jpg format

library(tidyverse)

figs <- "files/"  %>%  list.files(pattern = "tiff", full.names = T)

convert <- 1:length(figs) %>% map(\(x) {
  
  figname <- figs[x] %>% basename() %>% gsub("tiff", "jpg", .)
  
  figs[x] %>% 
    tiff::readTIFF(., native=TRUE) %>% 
    jpeg::writeJPEG(., target = file.path("files", figname))
  
})